Incendios en Villa Albertina - Sentinel 2
Visualización de imagenes de sentinel 2 respecto a los <a href='https://www.cadena3.com/noticia/radioinforme-3/continua-el-combate-del-incendio-en-villa-albertina_268066'>incendios</a> en Villa Albertina.
Data
Para determinar las imagenes de la zona afectada usamos los vectores provisto por la plataforma geojson.io.
Las imagenes satelitales las descargamos desde la plataforma de sentinel, via la api sentinelsat.
Las imágenes corresponden a los dias
- 20200815: _S2A_MSIL2A_20200815T141741_N0214_R010_T20JLM20200815T182929.zip
- 20200817: _S2B_MSIL2A_20200817T141049_N0214_R110_T20JLM20200817T163953.zip
- 20200820: _S2B_MSIL2A_20200820T141739_N0214_R010_T20JLM20200820T182855.zip
- 20200822: _S2A_MSIL2A_20200822T141051_N0214_R110_T20JLM20200822T182240.zip
# from http://geojson.io/#map=11/-30.5995/-64.2055
vector_file='data/villa_albertina.geojson'
# images from https://scihub.copernicus.eu using sentinelsat api
zip_20200815='data/S2A_MSIL2A_20200815T141741_N0214_R010_T20JLM_20200815T182929.zip'
zip_20200817='data/S2B_MSIL2A_20200817T141049_N0214_R110_T20JLM_20200817T163953.zip'
#
zip_20200820='data/S2B_MSIL2A_20200820T141739_N0214_R010_T20JLM_20200820T182855.zip'
zip_20200822='data/S2A_MSIL2A_20200822T141051_N0214_R110_T20JLM_20200822T182240.zip'
#collapse
class bands_etl(object):
"""
"""
@staticmethod
def get_ndi(band_1,band_2,weight=[1.0,1.0],dtype=rio.float32,nodata=-10000):
ndi = (weight[0]*band_1.astype(float)-weight[1]*band_2.astype(float))/(band_1*weight[0]*+band_2*weight[1])
ndi_val=np.where((band_1*weight[0]*+band_2*weight[1])>0,ndi,nodata)
return ndi_val.astype(dtype)
@staticmethod
def get_band_paths(path_,pattern_=['*_10m.jp2']):
"""
"""
names_bands={}
all_names_bands={}
for root, _, _ in os.walk(path_):
for p in pattern_:
ret=glob.glob(root+os.path.sep+p)
if len(ret)>0:
j=0
for r in ret:
full_path=os.path.join(root, r)
name=os.path.basename(full_path)
band=name.split('_')[-2]
c={band:full_path}
names_bands.update(c)
all_names_bands.update({p+'_'+str(j):full_path})
j+=1
else:
pass
return names_bands,all_names_bands
@staticmethod
def basic_resampling(input_raster,output_raster,resampling_method='cubic',cell_size=20):
"""
:param input_raster: file to be resampled
:param output_raster: file to be written
"""
warp_opts = [
"-r",
resampling_method,
"-tr",
str(cell_size),
str(cell_size),
]
#
gdal.UseExceptions()
ds = gdal.Warp(output_raster, input_raster, options=warp_opts) # noqa
del ds
def set_band_names(self,band_dict):
"""
:param band_dict: dict with band names and paths
"""
self.band_names_=band_dict
def get_band(self,b,open_flag):
"""
"""
if open_flag:
band_=rio.open(self.band_names_[b])
else:
band_=b
return band_
def set_mask_from_gpd(self,mask_gpd_geometry):
"""
"""
self.mask_gpd_geom_=mask_gpd_geometry
def get_clipping_params(self,band_init,mask_opts={'crop':True},open_flag=True):
"""
"""
band_=self.get_band(band_init,open_flag)
out_image_init, out_transform_init = mask(band_, self.mask_gpd_geom_,**mask_opts)
self.image_shape_0_=out_image_init[0].shape[0]
self.image_shape_1_=out_image_init[0].shape[1]
self.transform_=out_transform_init
self.dtype_=out_image_init[0].dtype
self.crs_=band_.crs
def bands_clipped_to_tiff(self,bands_list,tiff_name_='bands.tiff',driver_='Gtiff',open_flag=True,scale_band=False):
# Create an RGB image
if scale_band:
dtype_=rio.float32
else:
dtype_=self.dtype_
with rio.open(tiff_name_,'w',driver=driver_, width=self.image_shape_1_, height=self.image_shape_0_,
count=len(bands_list),crs=self.crs_,transform=self.transform_, dtype=dtype_) as img:
for i,b in enumerate(bands_list,start=1):
band_=self.get_band(b,open_flag)
out_image_init, out_transform_init = mask(band_, self.mask_gpd_geom_,crop=True)
if scale_band:
res=out_image_init[0]/out_image_init[0].max()
img.write(res.astype(dtype_),i)
else:
img.write(out_image_init[0],i)
#
band_.close()
out_image_init=None
out_transform_init=None
img.close()
def ndi_to_tiff(self,band1_,band2_,weight=[1.0,1.0],tiff_name_='ndi.tiff',driver_='Gtiff',open_flag=True,nodata=-10000):
# Create an NDVI image
with rio.open(tiff_name_,'w',driver=driver_, width=self.image_shape_0_, height=self.image_shape_1_,
count=1,crs=self.crs_,transform=self.transform_, dtype=rio.float32,nodata=nodata) as img:
band2_init=self.get_band(band2_,open_flag)
band1_init=self.get_band(band1_,open_flag)
band2, _ = mask(band2_init, self.mask_gpd_geom_,crop=True)
band1, _ = mask(band1_init, self.mask_gpd_geom_,crop=True)
band2_init.close()
band1_init.close()
img.write(self.get_ndi(band1,band2,weight,nodata=nodata))
img.close()
def write_array_single(
input_filename,
output_filename,
array,
array_dtype=rio.uint8,
nodata_=0,
attribs_source_="profile",
):
"""
:param input_filename: input filename (vrt,tiff,etc) to be taken as src for profile
:param output_filename: output filename
:param array: numpy array to be written
:param array_dtype: rasterio or numpy dtype
:param nodata_: nodata to be written as integer (e.g.:-10000)
:param attribs_source_: meta or profile
"""
with rio.open(input_filename) as src:
attribs_ = getattr(src, attribs_source_)
# update meta
attribs_.update(
dtype=array_dtype,
count=1,
driver="GTiff",
nodata=nodata_,
compresion="lzw",
)
with rio.open(output_filename, "w", **attribs_) as dst:
dst.write(array.astype(array_dtype), 1)
check_file = lambda path: os.remove(path) if os.path.isfile(path) else None
gpd_va=gpd.read_file(vector_file)
gpd_va_proj = gpd_va.to_crs(epsg=32720)
#
gpd_vt=gpd.read_file('data/villa_del_totoral.geojson')
gpd_vt_proj = gpd_vt.to_crs(epsg=32720)
#
cases_dict={}
#
file=zip_20200817
date_='20200817'
dirtmp=tempfile.mkdtemp()
with zipfile.ZipFile(file,"r") as zip_ref:
zip_ref.extractall(dirtmp)
#
# instanciamos el proc
bp=bands_etl()
# incluimos la capa vectorial
bp.set_mask_from_gpd(gpd_va_proj.geometry)
# Obtenemos los paths
paths,_=bp.get_band_paths(dirtmp, pattern_=['*B0?_10m.jp2','*B1?_20m.jp2'])
#
# Debido a las diferencias en el pixel size de la bandas debemos resamplear
resampled_bands={}
for band in paths:
case='_'+date_+'_Resampled'
name=os.path.join(products_tmp,band+case)
resampled_bands.update({band:name})
bp.basic_resampling(paths[band],name)
#
shutil.rmtree(dirtmp)
# seteamos las bandas
bp.set_band_names(resampled_bands)
# Generamos los parametros para el recorte
bp.get_clipping_params('B02')
#
case='false_color_urban_'+date_+'.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B12','B11','B04'],tiff_name_=tiff_name_,scale_band=True)
cases_dict.update({case:tiff_name_})
#
case='false_color_urban_'+date_+'_NOTScaled.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B12','B11','B04'],tiff_name_=tiff_name_,scale_band=False)
cases_dict.update({case:tiff_name_})
#
case='rgb_'+date_+'.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B04','B03','B02'],tiff_name_=tiff_name_,scale_band=True)
cases_dict.update({case:tiff_name_})
#
case='rgb_'+date_+'_NOTSCaled.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B04','B03','B02'],tiff_name_=tiff_name_,scale_band=False)
cases_dict.update({case:tiff_name_})
Ahora vamos a estimar (en forma muy aproximada) el area afectada. Para ello nos valemos de la banda 11. Vemos que es la que mayor contraste parece tener
Analizamos los valores y efectuamos el "thresholding" con el fin de quedarnos con una imagen en blanco y negro.
#collapse
file_path_threshold=os.path.join(products_tmp,'thresholding_b11_20200817.tif')
check_file(file_path_threshold)
!gdal_calc.py --calc="A>1400" -A 'products/false_color_urban_20200817_NOTScaled.tiff' --A_band=2 --outfile={file_path_threshold} --quiet
Vemos que si bien remarcamos la mayoria de la zona, quedan algunas extras. Las intentamos remover e invertimos la imagen.
#collapse
file_path_sieved=os.path.join(products_tmp,'thresholding_b11_sieved_20200817.tif')
check_file(file_path_sieved)
!gdal_sieve.py {file_path_threshold} {file_path_sieved} -st 10
Obtenemos el vector con la zona afectada via gdal-polygonize y visualizamos tanto el raster como el shapefile.
#collapse
vector_fire_sieved=os.path.join(products_tmp,'fire_shape_20200817.shp')
!gdal_polygonize.py {raster_to_shape} {vector_fire_sieved}
Estimamos nuevamente el area afectada
#
file=zip_20200822
date_='20200822'
dirtmp=tempfile.mkdtemp()
with zipfile.ZipFile(file,"r") as zip_ref:
zip_ref.extractall(dirtmp)
#
# instanciamos el proc
bp=bands_etl()
# incluimos la capa vectorial
bp.set_mask_from_gpd(gpd_va_proj.geometry)
# Obtenemos los paths
paths,_=bp.get_band_paths(dirtmp, pattern_=['*B0?_10m.jp2','*B1?_20m.jp2'])
#
# Debido a las diferencias en el pixel size de la bandas debemos resamplear
resampled_bands={}
for band in paths:
case='_'+date_+'_Resampled'
name=os.path.join(products_tmp,band+case)
resampled_bands.update({band:name})
bp.basic_resampling(paths[band],name)
#
shutil.rmtree(dirtmp)
# seteamos las bandas
bp.set_band_names(resampled_bands)
# Generamos los parametros para el recorte
bp.get_clipping_params('B02')
case='false_color_urban_'+date_+'.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B12','B11','B04'],tiff_name_=tiff_name_,scale_band=True)
cases_dict.update({case:tiff_name_})
#
case='false_color_urban_'+date_+'_NOTScaled.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B12','B11','B04'],tiff_name_=tiff_name_,scale_band=False)
cases_dict.update({case:tiff_name_})
#
case='rgb_'+date_+'.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B04','B03','B02'],tiff_name_=tiff_name_,scale_band=True)
cases_dict.update({case:tiff_name_})
#
case='rgb_'+date_+'_NOTSCaled.tiff'
tiff_name_=os.path.join(products,case)
bp.bands_clipped_to_tiff(['B04','B03','B02'],tiff_name_=tiff_name_,scale_band=False)
cases_dict.update({case:tiff_name_})
Observamos que todavia permanecen algunos focos activos.
Estimamos el area afectada
Alternativa para el computo de area afectada
En este caso podemos estimar el area via el NBR (Normalized Burn Ratio). Computamos el indice (_bands_etl().ndi_totiff). Efectuamos el thresholding poligonizamos y generamos el area ( en este caso para el dia 22/08)
